% $Id: e_pulse.m 437 2005-06-27 18:52:55Z edmaclin $
% e_pulse.m - a matlab function version of Gabriele's pulse correction algorithm
%
% 1) detect each beat (on a low-pass filtered copy of
%    the average (across channels) of the raw (un-normalized) DC data)
% 2) compute an average (across beats and channels) pulse waveform
%    for each data type (DC, PH, AC)
%    (requires 'warping' (interpolating) all beats to a common length)
% 3) 'detrend' the average waveform
% 4) compute regression of each beat (on each channel) with 'dewarped' average beat
%    to estimate contribution of pulse to total signal for each data type (DC, PH, AC)
% 5) subtract estimated pulse from raw data (amplitude=slope of regression)
%
% 6/13/02 - save average pulse for DC, AC & phase, & beat times
% 7/7/11 - read good_dist_chs, good_mrdc_s & good_mrac_chs from norm file
%          and only process good chs

function e_pulse(dood_dist_chs,log_file,plot_flg)

global hdr handles file_name i_block dc ac ph adc aac aph msg id_tags

warning off MATLAB:polyfit:RepeatedPointsOrRescale

tic

e_pulse_id = '$Id: e_pulse.m 437 2005-06-27 18:52:55Z edmaclin $';
check = strcmp(e_pulse_id,id_tags); % Make sure we don't already have the tag
if sum(check) == 0
    id_tags{size(id_tags,1)+1,1} = e_pulse_id;  % Add it in then.
end

% compute max and min beat period in sample points
%
MinBeat = round((hdr.sample_rate * 60) / hdr.max_BPM);
MaxBeat = round((hdr.sample_rate * 60) / hdr.min_BPM);

in_path=[hdr.data_path hdr.exp_name 'norm00-00\'];
out_path=[hdr.data_path hdr.exp_name 'pc00-00\'];

if ~strcmp(log_file,'')
    logid = fopen(log_file,'a');
else logid = 1; end

if ~exist(out_path)
    fprintf(logid,'Making apropriate directories...\n');
    dos(['mkdir ' out_path]); % !!! ??? won't compile...
end

if ~exist([out_path 'HR\'])
    dos(['mkdir ' out_path 'HR\']);
end

in_file_flg=fopen([in_path file_name]);
out_file_flg=fopen([out_path file_name]);

if in_file_flg > 0
	fclose(in_file_flg);
end
if out_file_flg > 0
	fclose(out_file_flg);
end

if ((in_file_flg>0) & (out_file_flg<0))
    
    scrsz = [1           1        1024         768];    %get(0,'ScreenSize');
    if plot_flg == 1
        fprintf(logid,'Plotting...\n');
        f3=get(0,'CurrentFigure');
    end
    
    %[status,msg] = mkdir(hdr.data_path, [hdr.data_type '00-00']);
    
    dc=0;ac=0;ph=0;mrdc=0;mrac=0;
    tmp_beat=0;
    
    % Load a data file
    %
    %file_name=[hdr.data_path hdr.data_type '00-00\' hdr.SS subj cf '.' block]
    fprintf(logid,'\ne-pulse reading %s\n',[in_path file_name]);
%     load([in_path file_name],'dc','ac','ph','mrdc','mrac');   
    load([in_path file_name],'-mat');
    
     if exist('good_mrac_chs')~=1
        % apply DC & AC quality thresholds
        %
        good_mrac_chs=find(mrac>=hdr.qual.norm_mrac_thr);
        good_mrdc_chs=find(mrdc>=hdr.qual.norm_mrdc_thr);
     end
    
    good_chs=intersect(good_mrac_chs,good_mrdc_chs);
    good_chs=intersect(good_chs,good_dist_chs);

    % zero the bad chans ( mostly for "old" normalized data)
    %
    temp=zeros(size(dc));
    temp(:,good_chs)=dc(:,good_chs);
    dc=temp;
    temp=zeros(size(dc));
    temp(:,good_chs)=ac(:,good_chs);
    ac=temp;
    temp=zeros(size(dc));
    temp(:,good_chs)=ph(:,good_chs);
    ph=temp;
    
    % 1) Detect those beats
    %
    %   First average the (un-normalized) DC data across channels
    %
    fprintf(logid,'\nAveraging raw data\n')
    AvgRawDC = zeros(hdr.n_points,1);                       % Initialize accumulator

    for i_ch=good_chs                               % For all good channels,
        RawDC(:,i_ch) = mrdc(i_ch).*dc(:,i_ch);       % Sum a de-normalized channel
    end

    AvgRawDC=mean(RawDC,2);
    AvgRawDC = AvgRawDC/length(good_chs);

    %   Filter the averaged DC data
    %
    %       Define a Butterworth filter
    %       Low & High pass freqs are specified as proportion of Nyquist freq
    %
    fprintf(logid,'Filtering averaged data\n')
    nf = hdr.sample_rate/2;                                 % Nyquist frequency
    filt_order=5; % use sharp cutoff (12) to remove breathing; 12 breaks filter!!!
    [b,a] = butter(filt_order,[hdr.BeatHiPass/nf,hdr.BeatLoPass/nf]);
    
    AvgFiltDC = filter(b,a,AvgRawDC);                   % Filter it
    
    %   NOW detect those beats
    %
    fprintf(logid,'Detecting beats: ')
    i_bc=1;                                             % Initialize beat counter
    CurMinBeat=1;                                       % Initialize first point of current window
    CurMaxBeat=MaxBeat-MinBeat;                 % Initialize last point of current window
    
    if plot_flg == 1
        f1=figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]); % Open a figure for plotting detected beats
    end

    while CurMaxBeat < hdr.n_points                         % Find all the beats
        
        [peak i_b]=max(AvgFiltDC(CurMinBeat:CurMaxBeat));       % Find peak (value & index) in current window
        Beat(i_bc)=CurMinBeat+i_b-2;                      % Record this peak's index in the "beat" array

        if plot_flg == 1
            hold on
            if i_bc+1 > 2
                offset=AvgFiltDC(Beat(i_bc-1));
                plot(AvgFiltDC(Beat(i_bc-1):Beat(i_bc))-offset);
                pause(.01);
            end
        end
        CurMinBeat=Beat(i_bc)+MinBeat;                  % Update first point of current window
        CurMaxBeat=Beat(i_bc)+MaxBeat;                  % Update last point of current window
        i_bc=i_bc+1;                                    % Increment beat counter
    end
    n_beats=i_bc-1;                                       % Record total number of beats
    
    fprintf(logid,'%4.0f beats detected!\n',n_beats);
    
    % 2) Compute the average warped pulse waveforms for each data type (DC, PH, AC)
    %
    fprintf(logid,'Warping and averaging\n')
    adc=zeros(hdr.n_warp,hdr.n_chans);                          % Initialize average pulse accumulators for DC,
    aph=zeros(hdr.n_warp,hdr.n_chans);                          %   Phase,
    aac=zeros(hdr.n_warp,hdr.n_chans);                          %   and AC
    
%     fprintf(logid,'ch) ');    
    for i_ch = good_chs                               % For each channel        
%         pause(.001);
%         fprintf(logid,'%3.0f ',i_ch);
%         if mod(i_ch,20) == 0
%             fprintf(logid,'\n');
%         end
        % DC first
        %
        for i_bc=4:n_beats                              % For each beat (Ignore incomplete 1st & last beats)
            clear tmp_beat                              % tmp_beat will vary in size
            tmp_beat=dc(Beat(i_bc-1):Beat(i_bc),i_ch);  % Stash current DC beat for interpolation
            
            % yi = interp1(x,y,xi,method)
            %   y is a vector containing the values of a function, and x is a vector of the same length
            %   containing the points for which the values in y are given. xi is a vector containing the
            %   points at which to interpolate. method is an optional string specifying an interpolation method
            %
            x=[1:length(tmp_beat)];
            xi=[1:hdr.n_warp];
            xi=(length(tmp_beat)-1)*xi/(hdr.n_warp)+1;
            fred = interp1(x,tmp_beat,xi,'linear');
            adc(:,i_ch)=adc(:,i_ch)+fred';
        end
        adc(:,i_ch)=adc(:,i_ch)/(n_beats-1);            % Normalize average beat (not actually necessary)
        % Phase
        %
        for i_bc=4:n_beats                              % For each beat (Ignore incomplete 1st & last beats)
            clear tmp_beat                              % tmp_beat will vary in size
            tmp_beat=ph(Beat(i_bc-1):Beat(i_bc),i_ch);  % Stash current Phase beat for interpolation
            x=[1:length(tmp_beat)];
            xi=[1:hdr.n_warp];
            xi=(length(tmp_beat)-1)*xi/(hdr.n_warp)+1;
            fred = interp1(x,tmp_beat,xi,'linear');
            aph(:,i_ch)=aph(:,i_ch)+fred';
        end
        aph(:,i_ch)=aph(:,i_ch)/(n_beats-1);            % Normalize average beat (not actually necessary)
        
        % AC
        %
        if i_ch == 256, fprintf('1 '); end
        for i_bc=4:n_beats                              % For each beat (Ignore incomplete 1st & last beats)
            clear tmp_beat                              % tmp_beat will vary in size
            tmp_beat=ac(Beat(i_bc-1):Beat(i_bc),i_ch);  % Stash current AC beat for interpolation
            x=[1:length(tmp_beat)];                     % x is just the index array
            xi=[1:hdr.n_warp];                              
            xi=(length(tmp_beat)-1)*xi/(hdr.n_warp)+1;      % xi are the points at which interpolation is desired
            fred = interp1(x,tmp_beat,xi,'linear');
            aac(:,i_ch)=aac(:,i_ch)+fred';
        end
        aac(:,i_ch)=aac(:,i_ch)/(n_beats-1);            % Normalize average beat (not actually necessary)
        
    end % i_ch

    % 3) 'detrend' the average waveforms (1st & last point = 0)
    %
%     fprintf(logid,'\nDetrending DC\nch) ')
    for i_ch=good_chs
%         pause(.001);       
%         fprintf(logid,'%3.0f ',i_ch);
%         if mod(i_ch,20) == 0
%             fprintf(logid,'\n');
%         end

        slope=(adc(hdr.n_warp-1,i_ch)-adc(1,i_ch))/(hdr.n_warp);    % Slope of line between 1st & last point
        for i=1:hdr.n_warp
            tmp(i,1)=adc(i,i_ch)-adc(1,i_ch)-slope*(i-1);   % Subtract trend line & offset (1st point)
        end
        adc(:,i_ch)=tmp;                                    % Recursion doesn't work here...
    end
    
    fprintf(logid,'\nDetrending AC\nch) ')
    for i_ch=good_chs
%         pause(.001);    
%         fprintf(logid,'%3.0f ',i_ch);
%         if mod(i_ch,20) == 0
%             fprintf(logid,'\n');
%         end
        
        slope=(aac(hdr.n_warp-1,i_ch)-aac(1,i_ch))/(hdr.n_warp);    % Slope of line between 1st & last point
        for i=1:hdr.n_warp
            tmp(i,1)=aac(i,i_ch)-aac(1,i_ch)-slope*(i-1);   % Subtract trend line & offset (1st point)
        end
        aac(:,i_ch)=tmp;                                    % Recursion doesn't work here...
    end
    
    fprintf(logid,'\nDetrending Phase\nch) ')
    for i_ch=good_chs
%         pause(.001);
%         fprintf(logid,'%3.0f ',i_ch);
%         if mod(i_ch,20) == 0
%             fprintf(logid,'\n');
%         end
        
        slope=(aph(hdr.n_warp-1,i_ch)-aph(1,i_ch))/(hdr.n_warp);    % Slope of line between 1st & last point
        for i=1:hdr.n_warp
            tmp(i,1)=aph(i,i_ch)-aph(1,i_ch)-slope*(i-1);   % Subtract trend line & offset (1st point)
        end
        aph(:,i_ch)=tmp;                                    % Recursion doesn't work here...
    end
    
    if plot_flg == 1
        close(f1);
        f1=figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]); % Open a figure for plotting avg pulses
        plot(adc);
        pause(.01);
    end
    
    % 4) compute regression of each beat (on each channel) with 'dewarped' average beat
    %    to estimate contribution of pulse to total signal for each data type (DC, PH, AC)
    %
    fprintf(logid,'\nComputing regression & subtracting\nch)')
    if plot_flg == 1
        f2=figure('Position',[scrsz(3)/2 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]); % Open a figure for plotting regression
    end
%     pcdc=zeros(hdr.n_points,hdr.n_chans);                        % Define 'pulse corrected DC' array
%     pcph=zeros(hdr.n_points,hdr.n_chans);                        % Define 'pulse corrected Phase' array
%     pcac=zeros(hdr.n_points,hdr.n_chans);                        % Define 'pulse corrected AC' array

    % copy raw data to create array for pulse corrected data (this leaves
    % "noise" in place of initial detected beats, instead of "0s" which
    % cause filtering problems
    %
    pcdc=dc;pcph=ph;pcac=ac;
    
    for i_ch=good_chs
%         pause(.001);        
%         fprintf(logid,'%3.0f ',i_ch);
%         if mod(i_ch,20) == 0
%             fprintf(logid,'\n');
%         end

        % DC
        %
        if plot_flg == 1
            hold on                        
        end
        for i_bc=4:n_beats
            clear tmp_beat xi x                         % tmp_beat & xi will vary in size
            tmp_beat=dc(Beat(i_bc-1):Beat(i_bc),i_ch);  % Stash current DC beat for interpolation
            xi=[1:length(tmp_beat)];
            xi=hdr.n_warp*xi/length(tmp_beat);              
            x=[1:hdr.n_warp];
            a_beat=interp1(x,adc,xi,'linear');          % Decimate average beat to size of current beat
            if sum(dc(:,i_ch)) ~= 0
                reg=polyfit(a_beat(:,i_ch),tmp_beat,1); % Compute regression
                if plot_flg == 1
                    plot(a_beat(:,i_ch),tmp_beat,'k.');      % Scatterplot avg vs. current beat
                    pause(.01);
                end
            else
                reg=[0 0];                              % For previously zeroed channels
            end
            
            % 5) subtract estimated pulse from raw data (amplitude=slope of regression)
            %
            pcdc(Beat(i_bc-1):Beat(i_bc),i_ch)=tmp_beat-reg(1)*a_beat(:,i_ch);
        end
        if plot_flg == 1
            hold off
        end
        % Phase
        %
        for i_bc=4:n_beats
            clear tmp_beat xi x                         % tmp_beat & xi will vary in size
            tmp_beat=ph(Beat(i_bc-1):Beat(i_bc),i_ch);    % stash current Phase beat for interpolation
            xi=[1:length(tmp_beat)];
            xi=hdr.n_warp*xi/length(tmp_beat);
            x=[1:hdr.n_warp];
            a_beat=interp1(x,aph,xi,'linear');
            if sum(dc(:,i_ch)) ~= 0
                reg=polyfit(a_beat(:,i_ch),tmp_beat,1);
                %plot(a_beat(:,i_ch),tmp_beat,'kd')
            else
                reg=[0 0];
            end
            
            % 5) subtract estimated pulse from raw data (amplitude=slope of regression)
            %
            pcph(Beat(i_bc-1):Beat(i_bc),i_ch)=tmp_beat-reg(1)*a_beat(:,i_ch);
        end
        
        % AC
        %
        for i_bc=4:n_beats
            clear tmp_beat xi x                         % tmp_beat & xi will vary in size
            tmp_beat=ac(Beat(i_bc-1):Beat(i_bc),i_ch);    % stash current AC beat for interpolation
            xi=[1:length(tmp_beat)];
            xi=hdr.n_warp*xi/length(tmp_beat);
            x=[1:hdr.n_warp];
            a_beat=interp1(x,aac,xi,'linear');
            if sum(dc(:,i_ch)) ~= 0
                reg=polyfit(a_beat(:,i_ch),tmp_beat,1);
                %plot(a_beat(:,i_ch),tmp_beat,'kd')
            else
                reg=[0 0];
            end
            
            % 5) subtract estimated pulse from raw data (amplitude=slope of regression)
            %
            pcac(Beat(i_bc-1):Beat(i_bc),i_ch)=tmp_beat-reg(1)*a_beat(:,i_ch);
        end
    end
    
    % Plot raw and pulse corrected data for comparison
    %
    if plot_flg == 1
        close(f3);                
        f3=figure('Position',[1 1 scrsz(3) scrsz(4)/2]); % Open a figure for plotting the results
        hold on
        for i_ch=good_chs
            
            plot(dc(:,i_ch)+.1*i_ch,'k');
            plot(pcdc(:,i_ch)+.1*i_ch+.05,'r');
            
            pause(.01);
        end
        xlim=([0 hdr.n_points]);
        hold off
        
        close(f1);
        close(f2);
        pause
    end
    
    % rename pulse corrected arrays to 'standard' names
    % (they will be distinguished by their directories)
    %
    dc=pcdc;
    ph=pcph;
    ac=pcac;
    
    % Write 'em out (if desired)
    %
    hdr_pulse=hdr; % save hdr WITH UNIQUE NAME !
    f_name=[out_path file_name];
    fprintf(logid,'Writing %s\n',f_name);
    if logid ~= 1 fprintf('Writing %s\n',f_name);end
    save(f_name,'dc','ph','ac','mrdc','mrac','id_tags',...
        'good_dist_chs','good_mrdc_chs','good_mrac_chs','hdr_pulse');

    % Save various heart rate data in "f_name_HR.nnn"
    %
    k=findstr('.',file_name);
    HR_f_name=[out_path 'HR\' file_name(1:k-1) '_HR' file_name(k:end)];
    save(HR_f_name,'AvgFiltDC', 'hdr', 'Beat','adc','aac','aph'); 

else
    if (in_file_flg<0)
        msg=sprintf('%s does not exist.\n',[in_path file_name]);
        fprintf(logid,'%s does not exist.\n',[in_path file_name]);
        fprintf('%s does not exist.\n',[in_path file_name]);
    end
    if (out_file_flg>0)
        msg=sprintf('%s already exists!\n',[out_path file_name]);
        fprintf(logid,'%s already exists!\n',[out_path file_name]);
        fprintf('%s already exists!\n',[out_path file_name]);
    end
end

if logid ~= 1, fclose(logid); end

toc

% $Log: e_pulse.m,v $
% Revision 1.32  2005/03/01 17:52:47  odellj
% Checking for pc00-00 directory was commented out due to compilation issues.  Uncommented since we're not compiling e_pulse at the moment.  Bug #27
%
% Revision 1.31  2005/01/24 17:21:01  odellj
% dc was doing the wrong channels
%
% Revision 1.29  2005/01/20 20:53:48  odellj
% fixed averaging to look only at the good_chs arrary, not every channel, hdr.n_chs.  Allows use of detector for LED pulse recording ala INM.
%
% Revision 1.28  2004/10/20 21:51:17  odellj
% Removed verbose statement
%
% Revision 1.27  2004/10/19 19:52:58  odellj
% Adjustments from INM -> PSC
%
% Revision 1.26  2004/09/09 21:06:09  odellj
% Made the statements to the screen friendlier when logging to file is checked.
%
% Revision 1.25  2004/09/07 19:10:26  odellj
% Cleaned up output.
%
% Revision 1.24  2004/09/07 16:31:23  edmaclin
% fixed fid tests
%
% Revision 1.22  2004/08/25 22:16:10  odellj
% Added the ability to pass the id_tag into the id_tags structure if this routine is called.  Also added some more checking in case directories were already made.  Using the exist command in regards to directories under DOS, which 'seems' to compile ok, unlike when used to check the existence of actual files.
%
% Revision 1.21  2004/08/25 18:38:09  odellj
% There was an errant pause statement outside of a plot loop.  Slipped it back in.
%
% Revision 1.20  2004/07/14 20:32:31  odellj
% changed the plot flag to be taken from the p_pod.fig
%
% Revision 1.19  2004/07/13 22:20:36  odellj
% Added the new variable _id to the subprogram to output to the inidiviually created files to identify the version of the producing function.
%
% Revision 1.18  2004/07/12 20:37:57  edmaclin
% *** empty log message ***
%
% Revision 1.17  2004/07/12 19:19:14  odellj
% typos, moved the opening of the log file closer to the top.
%
% Revision 1.16  2004/07/12 17:04:25  odellj
% Forgot to add the log_file parameter in the function declaration.
%
% Revision 1.15  2004/07/09 21:51:11  odellj
% Updated the logging of the subroutines.  Added log_file as a parameter for each of their calls.
%
% Revision 1.14  2004/07/07 18:08:32  odellj
% Prepped for logging.
%
% Revision 1.13  2004/07/01 18:08:50  edmaclin
% Fixed merged issues
%
% Revision 1.12  2004/06/30 18:14:19  odellj
% Added more verbose fprintf statements.
%
% Revision 1.11  2004/05/19 21:44:23  edmaclin
% good_chs not good_chans
%
% Revision 1.10  2004/05/19 20:36:23  edmaclin
% use .gch not .dis
%
% Revision 1.9  2004/05/18 15:44:02  edmaclin
% always writes output file (?? pc_save_flg=13??)
%
% Revision 1.8  2004/05/10 21:44:49  edmaclin
% declared SE arrays to avoid missing (n=0) bins
% process only "good" channels
%
% Revision 1.7  2004/05/07 22:16:16  edmaclin
% declared SD arrays to avoid missing (n=0) bins
%
% Revision 1.6  2004/04/06 19:16:20  edmaclin
% moved pause to plot loops
% vectorized mean beat computation
%
% Revision 1.5  2004/03/30 18:53:59  odellj
% Fixed the mkdir function.  Should compile now.
%
% Revision 1.4  2004/03/30 16:44:14  odellj
% Final(?) change of filt to norm.
%
% Revision 1.3  2004/03/29 22:11:33  odellj
% Changed norm to filt.
%
% Revision 1.1.1.1  2004/02/18 20:09:49  odellj
% initial import into CVS
%
% Revision 1.2  2004/02/05 20:49:24  odellj
% Added Id and Log tags.


